home *** CD-ROM | disk | FTP | other *** search
/ GFX Sensations 1 / Graphic Sensations - Volume 1.iso / tools / amiga / 3d_tools / irit40s.lha / Irit / cagd_lib / cbzr_aux.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-12-30  |  14.7 KB  |  467 lines
  1. /******************************************************************************
  2. * CBzr-Aux.c - Bezier curve auxilary routines.                      *
  3. *******************************************************************************
  4. * Written by Gershon Elber, Mar. 90.                          *
  5. ******************************************************************************/
  6.  
  7. #include <ctype.h>
  8. #include <stdio.h>
  9. #include <string.h>
  10. #include "cagd_loc.h"
  11.  
  12. /******************************************************************************
  13. * Given a bezier curve - subdivide it into two at the given parametric value. *
  14. * Returns pointer to first curve in a list of two curves (subdivided ones).   *
  15. * The subdivision is exact result of evaluating the curve at t using the      *
  16. * recursive algorithm - the left resulting points are the left curve, and the *
  17. * right resulting points are the right curve (left is below t).              *
  18. ******************************************************************************/
  19. CagdCrvStruct *BzrCrvSubdivAtParam(CagdCrvStruct *Crv, CagdRType t)
  20. {
  21.     CagdBType
  22.     IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
  23.     int i, j, l,
  24.     k = Crv -> Length,
  25.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
  26.     CagdRType
  27.     t1 = 1.0 - t;
  28.     CagdCrvStruct *LCrv, *RCrv;
  29.  
  30.     if (t < 0.0 || APX_EQ(t, 0.0) ||
  31.     t > 1.0 || APX_EQ(t, 1.0))
  32.     FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV);
  33.  
  34.     LCrv = BzrCrvNew(k, Crv -> PType);
  35.     RCrv = BzrCrvNew(k, Crv -> PType);
  36.  
  37.     /* Copy Curve into RCrv, so we can apply the recursive algo. to it.      */
  38.     for (i = 0; i < k; i++)
  39.     for (j = IsNotRational; j <= MaxCoord; j++)
  40.         RCrv -> Points[j][i] = Crv -> Points[j][i];
  41.  
  42.     for (j = IsNotRational; j <= MaxCoord; j++)
  43.     LCrv -> Points[j][0] = Crv -> Points[j][0];
  44.  
  45.     /* Apply the recursive algorithm to RCrv, and update LCrv with the         */
  46.     /* temporary results. Note we updated the first point of LCrv above.     */
  47.     for (i = 1; i < k; i++) {
  48.     for (l = 0; l < k - i; l++)
  49.         for (j = IsNotRational; j <= MaxCoord; j++)
  50.         RCrv -> Points[j][l] = RCrv -> Points[j][l] * t1 +
  51.                      RCrv -> Points[j][l+1] * t;
  52.     /* Copy temporary result to LCrv: */
  53.     for (j = IsNotRational; j <= MaxCoord; j++)
  54.         LCrv -> Points[j][i] = RCrv -> Points[j][0];
  55.     }
  56.     LCrv -> Pnext = RCrv;
  57.  
  58.     return LCrv;
  59. }
  60.  
  61. /******************************************************************************
  62. * Returns a new curve, identical to the original but with order N.          *
  63. *   Degree raise by multiplying by a constant 1 curve of order                  *
  64. * (NewOrder - Order).                                  *
  65. ******************************************************************************/
  66. CagdCrvStruct *BzrCrvDegreeRaiseN(CagdCrvStruct *Crv, int NewOrder)
  67. {
  68.     int i, j, RaisedOrder,
  69.     Order = Crv -> Order,
  70.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
  71.     CagdCrvStruct *RaisedCrv, *UnitCrv;
  72.  
  73.     if (NewOrder < Order) {
  74.     FATAL_ERROR(CAGD_ERR_WRONG_ORDER);
  75.     return NULL;
  76.     }
  77.     RaisedOrder = NewOrder - Order + 1;
  78.  
  79.     UnitCrv = BzrCrvNew(RaisedOrder, CAGD_MAKE_PT_TYPE(FALSE, MaxCoord));
  80.     for (i = 1; i <= MaxCoord; i++)
  81.     for (j = 0; j < RaisedOrder; j++)
  82.         UnitCrv -> Points[i][j] = 1.0;
  83.  
  84.     RaisedCrv = BzrCrvMult(Crv, UnitCrv);
  85.  
  86.     CagdCrvFree(UnitCrv);
  87.  
  88.     return RaisedCrv;
  89. }
  90.  
  91. /******************************************************************************
  92. * Return a new curve, identical to the original but with one degree higher    *
  93. * Let old control polygon be P(i), i = 0 to k-1, and Q(i) be new one then:    *
  94. *               i        k-i                          *
  95. * Q(0) = P(0), Q(i) = --- P(i-1) + (---) P(i), Q(k) = P(k-1).              *
  96. *               k         k                          *
  97. ******************************************************************************/
  98. CagdCrvStruct *BzrCrvDegreeRaise(CagdCrvStruct *Crv)
  99. {
  100.     CagdBType
  101.     IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
  102.     int i, j,
  103.     k = Crv -> Length,
  104.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
  105.     CagdCrvStruct
  106.     *RaisedCrv = BzrCrvNew(k + 1, Crv -> PType);
  107.  
  108.     for (j = IsNotRational; j <= MaxCoord; j++)                /* Q(0). */
  109.     RaisedCrv -> Points[j][0] = Crv -> Points[j][0];
  110.  
  111.     for (i = 1; i < k; i++)                        /* Q(i). */
  112.     for (j = IsNotRational; j <= MaxCoord; j++)
  113.         RaisedCrv -> Points[j][i] =
  114.         Crv -> Points[j][i-1] * (i / ((CagdRType) k)) +
  115.         Crv -> Points[j][i] * ((k - i) / ((CagdRType) k));
  116.  
  117.     for (j = IsNotRational; j <= MaxCoord; j++)                /* Q(k). */
  118.     RaisedCrv -> Points[j][k] = Crv -> Points[j][k-1];
  119.  
  120.     return RaisedCrv;
  121. }
  122.  
  123. /******************************************************************************
  124. * Return a normalized vector, equal to the tangent to Crv at parameter t.     *
  125. * Algorithm: pseudo subdivide Crv at t and using control point of subdivided  *
  126. * curve find the tangent as the difference of the 2 end points.              *
  127. ******************************************************************************/
  128. CagdVecStruct *BzrCrvTangent(CagdCrvStruct *Crv, CagdRType t)
  129. {
  130.     static CagdVecStruct P2;
  131.     CagdVecStruct P1, *T;
  132.     CagdBType
  133.     IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
  134.     int i, j, l,
  135.     k = Crv -> Length,
  136.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
  137.     CagdRType
  138.     t1 = 1.0 - t;
  139.     CagdCrvStruct *RCrv;
  140.  
  141.     if (APX_EQ(t, 0.0)) {
  142.     /* Use Crv starting tangent direction. */
  143.     CagdCoerceToE3(P1.Vec, Crv -> Points, 0, Crv -> PType);
  144.     CagdCoerceToE3(P2.Vec, Crv -> Points, 1, Crv -> PType);
  145.     }
  146.     else if (APX_EQ(t, 1.0)) {
  147.     /* Use Crv ending tangent direction. */
  148.     CagdCoerceToE3(P1.Vec, Crv -> Points, k - 2, Crv -> PType);
  149.     CagdCoerceToE3(P2.Vec, Crv -> Points, k - 1, Crv -> PType);
  150.     }
  151.     else {
  152.     if (t < 0.0 || t > 1.0)
  153.         FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV);
  154.  
  155.     RCrv = BzrCrvNew(k, Crv -> PType);
  156.  
  157.     /* Copy Crv into RCrv, so we can apply the recursive algo. to it. */
  158.     for (i = 0; i < k; i++)
  159.         for (j = IsNotRational; j <= MaxCoord; j++)
  160.         RCrv -> Points[j][i] = Crv -> Points[j][i];
  161.  
  162.     /* Apply the recursive algorithm to RCrv. */
  163.     for (i = 1; i < k; i++)
  164.         for (l = 0; l < k - i; l++)
  165.         for (j = IsNotRational; j <= MaxCoord; j++)
  166.             RCrv -> Points[j][l] = RCrv -> Points[j][l] * t1 +
  167.                        RCrv -> Points[j][l+1] * t;
  168.  
  169.     CagdCoerceToE3(P1.Vec, RCrv -> Points, 0, RCrv -> PType);
  170.     CagdCoerceToE3(P2.Vec, RCrv -> Points, 1, RCrv -> PType);
  171.  
  172.     CagdCrvFree(RCrv);
  173.     }
  174.  
  175.     CAGD_SUB_VECTOR(P2, P1);
  176.  
  177.     if (CAGD_LEN_VECTOR(P2) < SQR(EPSILON)) {
  178.     if (AttrGetIntAttrib(Crv -> Attr, "_tan") != TRUE) {
  179.         /* Try to move a little. This location has zero speed. However,  */
  180.         /* do it only once since we can be here forever. The "_tan"      */
  181.         /* attribute guarantee we will try to move EPSILON only once.    */
  182.         AttrSetIntAttrib(&Crv -> Attr, "_tan", TRUE);
  183.  
  184.         T = BzrCrvTangent(Crv, t < 0.5 ? t + EPSILON : t - EPSILON);
  185.  
  186.         AttrFreeOneAttribute(&Crv -> Attr, "_tan");
  187.  
  188.         return T;
  189.     }
  190.     else {
  191.         /* A zero length vector signals failure to compute tangent. */
  192.         return &P2;
  193.     }
  194.     }
  195.     else {
  196.     CAGD_NORMALIZE_VECTOR(P2);            /* Normalize the vector. */
  197.  
  198.     return &P2;
  199.     }
  200. }
  201.  
  202. /******************************************************************************
  203. * Return a normalized vector, equal to the binormal to Crv at parameter t.    *
  204. * Algorithm: pseudo subdivide Crv at t and using 3 consecutive control point  *
  205. * (p1, p2, p3) of subdivided curve find the bi-normal as the cross product of *
  206. * the difference (p1 - p2) and (p2 - p3).                      *
  207. *   Since a curve may have not BiNormal at inflection points or if the 3      *
  208. * points are colinear, NULL will be returned at such cases.              *
  209. ******************************************************************************/
  210. CagdVecStruct *BzrCrvBiNormal(CagdCrvStruct *Crv, CagdRType t)
  211. {
  212.     static CagdVecStruct P3;
  213.     CagdVecStruct P1, P2;
  214.     CagdBType
  215.     IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
  216.     int i, j, l,
  217.     k = Crv -> Length,
  218.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
  219.     CagdRType
  220.     t1 = 1.0 - t;
  221.     CagdCrvStruct *RCrv;
  222.  
  223.     /* Can not compute for linear curves. */
  224.     if (k <= 2)
  225.     return NULL;
  226.  
  227.     /* For planar curves, B is trivially the Z axis. */
  228.     if (CAGD_NUM_OF_PT_COORD(Crv -> PType) == 2) {
  229.     P3.Vec[0] = P3.Vec[1] = 0.0;
  230.     P3.Vec[2] = 1.0;
  231.     return &P3;
  232.     }
  233.  
  234.     if (APX_EQ(t, 0.0)) {
  235.     /* Use Crv starting tangent direction. */
  236.     CagdCoerceToE3(P1.Vec, Crv -> Points, 0, Crv -> PType);
  237.     CagdCoerceToE3(P2.Vec, Crv -> Points, 1, Crv -> PType);
  238.     CagdCoerceToE3(P3.Vec, Crv -> Points, 2, Crv -> PType);
  239.     }
  240.     else if (APX_EQ(t, 1.0)) {
  241.     /* Use Crv ending tangent direction. */
  242.     CagdCoerceToE3(P1.Vec, Crv -> Points, k - 3, Crv -> PType);
  243.     CagdCoerceToE3(P2.Vec, Crv -> Points, k - 2, Crv -> PType);
  244.     CagdCoerceToE3(P3.Vec, Crv -> Points, k - 1, Crv -> PType);
  245.     }
  246.     else {
  247.     if (t < 0.0 || t > 1.0)
  248.         FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV);
  249.  
  250.     RCrv = BzrCrvNew(k, Crv -> PType);
  251.  
  252.     /* Copy Crv into RCrv, so we can apply the recursive algo. to it. */
  253.     for (i = 0; i < k; i++)
  254.         for (j = IsNotRational; j <= MaxCoord; j++)
  255.         RCrv -> Points[j][i] = Crv -> Points[j][i];
  256.  
  257.     /* Apply the recursive algorithm to RCrv. */
  258.     for (i = 1; i < k; i++)
  259.         for (l = 0; l < k - i; l++)
  260.         for (j = IsNotRational; j <= MaxCoord; j++)
  261.             RCrv -> Points[j][l] = RCrv -> Points[j][l] * t1 +
  262.                        RCrv -> Points[j][l+1] * t;
  263.  
  264.     CagdCoerceToE3(P1.Vec, RCrv -> Points, 0, RCrv -> PType);
  265.     CagdCoerceToE3(P2.Vec, RCrv -> Points, 1, RCrv -> PType);
  266.     CagdCoerceToE3(P3.Vec, RCrv -> Points, 2, RCrv -> PType);
  267.  
  268.     CagdCrvFree(RCrv);
  269.     }
  270.  
  271.     CAGD_SUB_VECTOR(P1, P2);
  272.     CAGD_SUB_VECTOR(P2, P3);
  273.  
  274.     CROSS_PROD(P3.Vec, P1.Vec, P2.Vec);
  275.  
  276.     if ((t = CAGD_LEN_VECTOR(P3)) < EPSILON)
  277.     return NULL;
  278.     else
  279.     CAGD_DIV_VECTOR(P3, t);                /* Normalize the vector. */
  280.  
  281.     return &P3;
  282. }
  283.  
  284. /******************************************************************************
  285. * Return a normalized vector, equal to the normal to Crv at parameter t.      *
  286. * Algorithm: returns the cross product of the curve tangent and binormal.     *
  287. ******************************************************************************/
  288. CagdVecStruct *BzrCrvNormal(CagdCrvStruct *Crv, CagdRType t)
  289. {
  290.     static CagdVecStruct N, *T, *B;
  291.  
  292.     T = BzrCrvTangent(Crv, t);
  293.     B = BzrCrvBiNormal(Crv, t);
  294.  
  295.     if (T == NULL || B == NULL)
  296.     return NULL;
  297.  
  298.     CROSS_PROD(N.Vec, T -> Vec, B -> Vec);
  299.  
  300.     CAGD_NORMALIZE_VECTOR(N);                /* Normalize the vector. */
  301.  
  302.     return &N;
  303. }
  304.  
  305. /******************************************************************************
  306. * Return a new curve, equal to the derived curve.                  *
  307. * Let old control polygon be P(i), i = 0 to k-1, and Q(i) be new one then:    *
  308. * Q(i) = (k - 1) * P(i+1) - P(i), i = 0 to k-2.                      *
  309. ******************************************************************************/
  310. CagdCrvStruct *BzrCrvDerive(CagdCrvStruct *Crv)
  311. {
  312.     CagdBType
  313.     IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
  314.     int i, j,
  315.     k = Crv -> Length,
  316.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
  317.     CagdCrvStruct *DerivedCrv;
  318.  
  319.     if (!IsNotRational)
  320.     return BzrCrvDeriveRational(Crv);
  321.  
  322.     DerivedCrv = BzrCrvNew(MAX(1, k - 1), Crv -> PType);
  323.  
  324.     if (k >= 2) {
  325.     for (i = 0; i < k - 1; i++)
  326.         for (j = IsNotRational; j <= MaxCoord; j++)
  327.         DerivedCrv -> Points[j][i] =
  328.             (k - 1) * (Crv -> Points[j][i+1] - Crv -> Points[j][i]);
  329.     }
  330.     else {
  331.     for (j = IsNotRational; j <= MaxCoord; j++)
  332.         DerivedCrv -> Points[j][0] = 0.0;
  333.     }
  334.  
  335.     return DerivedCrv;
  336. }
  337.  
  338. /******************************************************************************
  339. * Return a new Bezier curve, equal to the integral of the given Bezier curve. *
  340. * The given Bezier curve should be nonrational.                      *
  341. *                                          *
  342. *           n           n           n       n+1                  *
  343. *   /         /-           -      /       -   P    -                  *
  344. *  |        | \        n       \     |  n       \    i   \  n+1              *
  345. *  | C(t) = | / P  B (t) = / P   | B (t) = / -----  / B   (t) =              *
  346. * /        /  -     i  i       -  i /   i       - n + 1  -  j              *
  347. *            i=0          i=0             i=0     j=i+1                  *
  348. *                                          *
  349. *        n+1 j-1                                  *
  350. *         -   -                                      *
  351. *     1   \   \        n+1                                  *
  352. * = ----- /   / P  B   (t)                              *
  353. *   n + 1 -   -  i  j                                  *
  354. *        j=1 i=0                                  *
  355. *                                          *
  356. ******************************************************************************/
  357. CagdCrvStruct *BzrCrvIntegrate(CagdCrvStruct *Crv)
  358. {
  359.     int i, j, k,
  360.     n = Crv -> Length,
  361.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
  362.     CagdCrvStruct *IntCrv;
  363.  
  364.     if (CAGD_IS_RATIONAL_CRV(Crv))
  365.     FATAL_ERROR(CAGD_ERR_RATIONAL_NO_SUPPORT);
  366.  
  367.     IntCrv = BzrCrvNew(n + 1, Crv -> PType);
  368.  
  369.     for (k = 1; k <= MaxCoord; k++) {
  370.     CagdRType
  371.         *Points = Crv -> Points[k],
  372.         *IntPoints = IntCrv -> Points[k];
  373.  
  374.     for (j = 0; j < n + 1; j++) {
  375.         IntPoints[j] = 0.0;
  376.         for (i = 0; i < j; i++)
  377.             IntPoints[j] += Points[i];
  378.         IntPoints[j] /= n;
  379.     }
  380.     }
  381.  
  382.     return IntCrv;
  383. }
  384.  
  385. /******************************************************************************
  386. * Convert a Bezier curve into Bspline curve by adding an open knot vector.    *
  387. ******************************************************************************/
  388. CagdCrvStruct *CnvrtBezier2BsplineCrv(CagdCrvStruct *Crv)
  389. {
  390.     CagdCrvStruct *BspCrv;
  391.  
  392.     if (Crv -> GType != CAGD_CBEZIER_TYPE) {
  393.     FATAL_ERROR(CAGD_ERR_WRONG_CRV);
  394.     return NULL;
  395.     }
  396.  
  397.     BspCrv = CagdCrvCopy(Crv);
  398.  
  399.     BspCrv -> Order = BspCrv -> Length;
  400.     BspCrv -> KnotVector = BspKnotUniformOpen(BspCrv -> Length,
  401.                                BspCrv -> Order, NULL);
  402.     BspCrv -> GType = CAGD_CBSPLINE_TYPE;
  403.     return BspCrv;
  404. }
  405.  
  406. /******************************************************************************
  407. * Convert a Bspline curve into a set of Bezier curves by subdividing the      *
  408. * Bspline curve at all its internal knots.                      *
  409. *   Returned is a list of Bezier curves.                      *
  410. ******************************************************************************/
  411. CagdCrvStruct *CnvrtBspline2BezierCrv(CagdCrvStruct *Crv)
  412. {
  413.     int i,
  414.     Order = Crv -> Order,
  415.     Length = Crv -> Length;
  416.     CagdRType LastT,
  417.     *KnotVector = Crv -> KnotVector;
  418.     CagdCrvStruct
  419.     *BezierCrvs = NULL,
  420.     *OrigCrv = Crv;
  421.  
  422.     if (Crv -> GType != CAGD_CBSPLINE_TYPE) {
  423.     FATAL_ERROR(CAGD_ERR_WRONG_CRV);
  424.     return NULL;
  425.     }
  426.  
  427.     for (i = Length - 1, LastT = KnotVector[Length]; i >= Order; i--) {
  428.         CagdRType
  429.             t = KnotVector[i];
  430.             
  431.     if (!APX_EQ(LastT, t)) {
  432.             CagdCrvStruct
  433.         *Crvs = BspCrvSubdivAtParam(Crv, t);
  434.  
  435.             if (Crv != OrigCrv)
  436.                 CagdCrvFree(Crv);
  437.  
  438.             Crvs -> Pnext -> Pnext = BezierCrvs;
  439.             BezierCrvs = Crvs -> Pnext;
  440.  
  441.             Crv = Crvs;
  442.             Crv -> Pnext = NULL;
  443.  
  444.         LastT = t;
  445.         }
  446.     }
  447.  
  448.     if (Crv == OrigCrv) {
  449.     /* No interior knots in this curve - just copy it: */
  450.     BezierCrvs = CagdCrvCopy(Crv);
  451.     }
  452.     else {
  453.         Crv -> Pnext = BezierCrvs;
  454.         BezierCrvs = Crv;
  455.     }
  456.  
  457.     for (Crv = BezierCrvs; Crv != NULL; Crv = Crv -> Pnext) {
  458.         Crv -> GType = CAGD_CBEZIER_TYPE;
  459.     Crv -> Length = Crv -> Order;
  460.     IritFree((VoidPtr) Crv -> KnotVector);
  461.     Crv -> KnotVector = NULL;
  462.     }
  463.  
  464.     return BezierCrvs;
  465. }
  466.  
  467.